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Abstract: A four-parameter quadrilateral nonconforming finite element with DSP (double set pa- 
rameters) is proposed to approximate the pure displacement linear elasticity problem for 
homogeneous isotropic material. It is shown that the performance of the scheme does not 
deteriorate as the material becomes nearly incompressible. Numerical results are also given 
to verify the theoretical analysis. Furthermore, the superconvergence result is obtained 
through a postprocessing operator. All the results can be extended to three dimensional 
case. 
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1 Introduction 


There now exists a large number of publications devoted to the task of constructing and 
analyzing finite element approximations for problems in solid mechanics, in which it is necessary 
to circumvent volumetric locking. Interested readers can refer to [1-4]. 

In order to avoid the locking effects of the nearly incompressible linear elasticity, a new 
nonconforming scheme is presented in this paper. And the nonconforming finite element (named 
as QB element) is constructed by the DSP (double set parameters) method. Then we investigate 
the quadrilateral QB nonconforming scheme for the linear elasticity problem. For rectangular 
meshes, we can equate the nonconforming method with a mixed method directly. As to the 
equivalent mixed formulation, we find that it exhibits some features of the well-known Q- 
element for the Stokes problem. Following the ideas presented in [5-8], and by using some special 
relations between QB element and Q; element, the Babuška-Brezzi conditionl®! with a constant 
(inf-sup constant) independent of h is derived after filtering out some “local spurious” pressure 
modes. Then an optimal convergence rate is established for both displacement and smoothed 
pressure. Based on the mathematical analysis of the QB — Po scheme for the mixed formulation, 
and by introducing another modification of the mixed form as an auxiliary problem, the optimal 
error estimate in energy norm is obtained. Finally, using the interpolation postprocessing 
technique, the superconvergence relult is obtained. 
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2 Notations and preliminaries 


In the context of elasticity, vector- and tensor- or matrix-valued functions are written in 
boldface form. Let Q C R? be a convex polygonal domain. We use standard definitions for the 
Sobolev spaces H*(Q) and their associated norms ||- ||5,0, and seminorms |- |s o, and 


LEQ) = fa E IO); | ga = o}. 


Finally, P,(Q) denotes the space of all polynomials on Q of degree < k, and Q}(N) denotes the 
space of all polynomials on Q of degree < k in each variable. 
The pure displacement planar elasticity problem for homogeneous isotropic material can be 


described as 
—pAu — (u + à)graddivu =f, in Q, 


(1) 
u= 0, on OQ, 
where u denotes the displacements, f the body forces, u and A are the Lamé parameters. 
The corresponding variational problem is 
Find u c H}(Q), such that (2) 
2 


u(7u, Vv) + (A + 2) (divu, divv) = (f, v), Vv € Hj(Q). 


The following theorem yields solvability of problem (2), which can be found in [2,9]. 
Theorem 2.1 Let Q C R? be a bounded convex polygonal domain, f € L?(Q), then the 
variational form (2) has a unique solution u € Hj(Q) N H?(Q), and there exists a positive 


constant C such that 
[ull2,.0 + Alldivullio < Cllfllo,o- (3) 


If we define p = (A+ y:)divu, then the equations of elasticity may also be written in a mixed 
form: find (u, p) € Hj() x £3(Q) such that 


u(Vu, Vv) + (divv, p) = (fv), Vv € Hi), 
(A+ u)! (p,q) — (divu,g) =0, Yg € L2(Q). 


(4) 


This formulation is valid even in the incompressible limit A = oo. For the application in fluid 
flow, ys is the viscosity, and u is the velocity. At the incompressible material limit A = 00, 
we obtain from problem (4) the equation for incompressible elasticity, the same as the Stokes 
equation. 


3 Double set parameter method and QB element 


As is known that the main disadvantage of nonconforming elements is more total degrees of 
freedom involved than conforming ones. Taking for the rotated Q; element!) as an example, 
its total degree of freedom is almost double as that of conforming bilinear element. -A distinct 
idea is to combine the two elements together, and then the double set parameter method! 
follows it. The main advantage of this new approach is the nodal parameters and the degree 


120 CHINESE JOURNAL OF ENGINEERING MATHEMATICS VOL. 28 


of freedom may be chosen independently. In principle, the nodal parameters are chosen to be 
simple so that the total number of unknowns in the resulting discrete system is small, while 
the degrees of freedom are selected to meet convergence requirement setting by the generalized 
patch test!!). Then, following [12], we will present a nonconforming element, which is named 


QB element. 
Assume that K = [—1, 1] x [—1, 1] is the reference element with vertices 


@=(-1,~1), @=(1,-1), @ = (1,1), & = (—1,1), 


and sides 


~ 


1=@@2, I2=G2G3, l3 =G3G4, l4 = G4ai. 

The first set of degrees of freedom on K is taken as the mean-value-oriented rotated Qı element 
i ect Ae eel a NE 

D(@) = (012, 593, 034, da1) , (5) 


where i 
Vii+1) = at ds, i= 1, 2,3, 4. 
i i 


The shape function space on reference K is taken as 
P = P\(K)U {6 - 17}. (6) 


Suppose 
T = By + B£ + Ban + Palé? — n’). (7) 


Substitute (7) into (5) one can get 


bı = 402 + 023 +334 + G41), Bo = 3 (23 — da), 


(8) 
b3 = 3 (G34 —%12), Ba = 3(-T12 + Dog — D34 + a1). 
Using the double set parameters method, we take another set of nodal parameters 
A Ana a AAT 
QO) = (G1, d2, 33, d4) , (9) 
i.e., the function values at 4 nodes, which are the real degrees of freedom. 
Approximating the degrees of freedom D(F) by the trapezoidal rule as follows 
p ae Pes A 
Vitis) = gv + ¥i41)+46:(@), 4 =1,2, 3,4. (10) 


For any v € Py (R ), these discretizations are exact, and an application of the Bramble-Hilbert 
lemma gives 6;(0) < CJ]; g- 
Suppose that the affine mapping Fx : K — K is defined as 
z= TK + hgg, 


Y = yk + hgn. 


NO. 1 Sun Huixia, Xiao Liuchao: Locking-free DSP Element for Linear Elasticity Problem 121 








On the physical element K, the above process can be expressed as 
D(v) = GQ(v) + d(v), (11) 
where 


Div) = Do Fx"), Q(v)=QGoFg'), dv) = (61(v), 62(v), 63(v), da(v)) *, 


3 4 0 0 

a =i ; 0 3 3 0 
6;(v) = 6:(®@)o Fx’? =O(h)\vlaxn, i=1,2,3,4, G= 

0033 

200 4 


Dropping 5; (8) at the right hand of (10), then we get the expressions of 8; of (7) in Q(%): 


(12) 


| bi =} ++ +0), b2 = i0 +D -D — hy), 
Bs = (03 +04 -D — 32), pa =0. 
The shape function space on K is defined by 
Pg = {v=o Fz}, 3e P is determined by (7) and (12)}. 
Assume T and I x to be the finite element interpolation operators deduced by P and Pr, 
respectively, then Inv = (T 0) o Fg’, where 
To = pı + b£ + Ban, Vie H?(K), (13) 
with the forms of (1, G2, @3 defined as (12). 
The finite element space Vp associated with QB element is given by 
Vn = {un; valk =D Fj’, © is defined by (13), VK € Jh, 
v(a) =0, for any vertex a on ôN}, (14) 


where Jn is a subdivision of Q. 

Let In : C°(Q) — Vn be the global interpolation of QB element on 2 with In| = Ix. 

Remark 3.1 It follows from (13) that the interpolation of QB element is the same as 
the linear part of conforming bilinear element. In fact, one can check that the interpolation of 
QB element is a disturbance of that of bilinear conforming element, i.e., (Ixv)(ai) = v(ai) + 
O(hx)|vl2,x- So it is a nonconforming element, but this does not influence the good properties 
of QB element as follows. 

Lemma 3.1 (i) Let Va be the bilinear finite element space, i.e., 


Vin = {v € ©°(Q); vle E€ QK), VK € Jn, vlon =O}, 
then we have 
Va = Vin = {7, Jo Fg is the linear part of vo Fg, VK En, vE Vin}. (15) 
(ii) For the interpolation operator I», there holds the error estimate 


lv — Invilo.a + Ally —Inulln < Ch llulian, Vv € H?N). (16) 
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4 Error estimates 


In this section we consider the nonconforming QB scheme based on the variational formu- 
lation (2). 

Let Vn = Vn x Vh, Von = Von X Von- For vn € Vn, we define Viva to be the L?() function 
whose restriction to each quadrilateral K € J, is given by Vva|g. Analogous definitions hold 
for div, and rotp. 

Define the discrete semi-norm 


1 
2 
Ivallma = (SD valar ) 


KETn 


The nonconforming finite element approximation scheme for problem (2) is then given as follows 
Find u, E Va, such that 
(17) 
L(VhUh, Viva) + (A+ u)(divrun, divava) = (f, va), Yva E Va. 


Remark 4.1 Brenner!#! analyzed the nonconforming linear triangular finite element meth- 
ods for the problem (1). In the nonconforming case, they used the mid-point values on each 
edge of elements as degrees of freedom, and thus as compared with the finite elements here, the 
degrees of freedom needed in [2] are about three times of ours. 

Remark 4.2 Recently, [4] proposed a nonconforming rectangular element following the 
ideas of Brenner]. It is a Pz rectangular element with the interpolation operator preserving 
zero-divergence at each element. Because of the large degrees of freedom, complex construction 
of this element and nonapplication of more general meshes, it is not an advisable locking free 
element. 

Now we introduce a piecewise constant finite element space 


Mn = {an € r3 (9); anlk € Po(K), VK € In}. 


As a counterpart of the problem (17), we introduce the following mixed formulation: find 
(un, ph) E Va X My, such that 
L(JhUh, Vava) + (divava, pr) = (E, Va), Vvan € Va, (18) 
(A+ u) (pai dn) — (divaun, gn) = 0, Van E€ Mn. 
Then it can be easily proved that the finite element approximation problem (17) is equivalent 
to the mixed problem (18) with 


Palk = (At+p)divaunlk, VK € Tn. 


Let us stress that the introduction of the new variable is done only as a tool for the math- 
ematical analysis of the method implemented by (17). Based on the above discussion, we can 
concentrate on the mixed problem (18). In order to establish the uniform convergence with 
respect to À, it is sufficient to consider the limiting case when 4 —> oo, which leads to 

H(Wrun, VhVh) + (divavn, ph) = (f, va), Vva € Va, (19) 
(divnun, gh) =0, Van E Mn. 
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For the sake of convenience, we assume that the subdivision J, has been obtained from an 
arbitrary regular affine quadrilateral meshes Jo, = {M} by dividing each element of Jan into 
four congruent rectangles. It can be easily seen that the pair (Vr, My) does not always satisfy 
the inf-sup condition. In fact, for any v € Vpn, we have V € Vp, and 


f qndivvdady = J qndivvdzdy, if qndivVdady = f qndivvdady, Y qnr E€ My. 
K K 2 re) 


So the pair (Vn, Mp) shares the same stable properties of (Ven, Mn) under affine quadrilateral 
meshes. It is well-known that the checkerboard function may cause (Von, Mh) a spurious 
pressure mode on a uniform grid. Interested readers can refer to [5-7] for detailed analysis. 

We are working with a finite dimensional space Ph © Mp, such that the pair (Va, Ph) 
satisfies the inf-sup condition. In order to exclude the possibility of spurious or nearly spurious 
pressure modes, the local basis functions for P, on a 2 x 2 patch of M are taken as i,m, which 
are indicated in Figure 1. Thus we include in P, all piecewise constants of the form 


3 
P, = {an E LS(Q) : qalm = XJ aim Yim, Me Jan}. 
i=l 





























P1,M P2,M P3,M 
Figure 1: Local basis functions of Ph 


Next, the inf-sup condition for the mixed problem (19) is established. 
Lemma 4.1 Under the above hypothesis, the pair (V;,, Pp) satisfies the following Babuška- 
Brezzi condition or inf-sup condition! 


b 
sup AON) > Cation, Van € Px (20) 


OfvnEV), Ivan llr, 


where C is a positive constant independent of h. 
Proof By the references [5-7], the pair (Ven, Pp) satisfies a uniform inf-sup condition, i.e., 


sup bh (v, dh) 


2C\lanlloe, Van € Ph- (21) 
OAVEVon lvii, 


An equivalent formulation to (21) is: for any gn € Ph, there is a function v € Vpn such that 


br{v, dn) > Cilgallo livla < Cllanllo.a. (22) 
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Then taking V € Vn, by the scaling argument and inverse inequality, one can conduct 


lv — Vin < Chiiviian < Clivina: (23) 

Thus we have obtained 
bn(¥, qn) = ba (v, gn) = Cllanllo,a, (24) 
I¥llan < Cllvllia < Cllanllo,a, (25) 


which implies the desired assertion. 
Let ph € P, be the L? projection of pa on Ph. Then the pair (ua, Ph) € (Va, Pr) is the 
unique solution of the following equations 


H(Vrun; Vava) + (divava, Ph) = (f, Va), V Vhr € Va, a 
(divaun, qn) = 0, Vv dn E Pp. 


Now, we are in the position to derive the error estimates. 
Theorem 4.1 Let (u, p) be the solution of problem Mo» or (4) and (up, pp) be the solution 
of problem (18) or (19). Then we have 


lu — uallija + lip — Pallo,e < CAllfllo,a. (27) 
Proof Owing to [6], the following error estimate holds 
llu — unlli,r + Ilp — Pallo,a 
< c{ inf |u- v inf ||p — Gp 
< Cy inf Il allan + ui lp — fallon 


+ sup u(Jhu, VaVn) as (diva Vh, p) z (£, Vh) }. (28) 
vnéVa IlValli.n 


We only need to bound the three terms at the right hand of (28) one by one. 
Suppose M = Uni Ki € Jan with K; € Jn (1 < i < 4). For the sake of the subsequent 
analysis, the operator jn is defined by 


. Pok:p + łam, i = 1,3, 
JhP = (29) 
Pox,p = dam, i = 2,4, 


where am = Pox,p+ Pox,p — Pok, p — Poxsp- 
A direct calculation shows that japlm = &1,MP1,M + &2,MP2,M + &3,MP3,M, Where 


1 

Q1,M = ekp + Pokap + Poxsp + Pok,P), 
1 

Q2, M = qPox.P — Pox,p — Pox,p + Pox,p); 


1 
a3,M = q(-Por.P ~ Pox,p + Pox3p + Por,p). 
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This implies that for any p € L2(Q), jap € Ph- 
Let Qj7 = am © Fx, it can be checked easily that 


ân < CllPlo m <CllAl.g@, ĉm=0 YPE Po(M). (30) 
An application of Bramble-Hilbert Lemma yjeldsl!3] 
lom| < CP a < Clph,m, Vpe H'(M). (31) 


Then there holds 
i 


è ~ : : 2 
inf lp- ilon < lp- jploa=( Y lp- jnvl8,ne ) 
gan EPh MEJan 


4 i 
( E So (lp- Porplê x + helen?) )* <Chlphia. (32) 


METan t=1 


IA 


Let I» = In x In, and by (16), there holds 


ying, Iu- valia < lu — Trulla < Chlfulla,a- (33) 


f vras= [vas = f vas 
F F F 


is continuous along edges of elements and vanishes if F C 00, where v € Van, F C OK, for all 
K € Jn, using the technique developed by [2,9], the consistency error can be estimated as 


Since 


L(VhU, Vavn) + (divava, p) — (£, va) < Ch(lulle,0 + lela) Il vallsn- (34) 


A combination of (32)-(34) and the regularity result (3) follows (27), which completes the proof. 

In the subsequent part of this section we will show that the nonconforming finite element 
solution up, defined by (17) converges to the solution u satisfying (2) in energy norm. The 
energy norm is defined by 





Il lla = Ve(Van-s Vr) + (A + 2)(diva:, diva) . 


Now, we will turn back to the primitive displacement variable approximation problem (17). 
The famous second Strang’s Lemma!:!5) gives 


; |En(u, Va)| 
u—u <C inf |u- vala + su =}, 35 
lu- usla <C{ inf, u- valat sa oT} (35) 


where 
E,,(u, vah) = L(JhU, Travan) + (A+ »)(divnu, divava) — (f, va). 


Similar to (34), the second term at the right hand of (35), i.e., the consistency error can be 
estimated as 


|En(u, va)| < Ch([lullao + (A + p)lldivulh, a) llvali a < Callflloallvalln. (36) 
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Let us consider a slight modification of (18): find (ŭn, Da) € Vn x My such that 


H(VnUn, Vavn) + (divava, Bn) = (f, Va), Viva € Va; 


(37) 
(A + Bw) (p, gn) koz (divrùn, gn) = 0, v qh E€ Mp- 
Based on the mathematical analysis of problem (19) and (18), one can obtain that 
ju — Ballin < Chllfllo,o, (38) 
div, ty, = (A+ u) + Ponp = Po,divu, (39) 
where Pon is the L? projection operator on Mp. 
Then by taking va = tt, in (35), we have 
lu- tall, = xllu—dall?, + A+ p)lidivu — divrùz ll o 
< Ch*|[fll,o + Ch?(A + p)||divull7.g < CRIEI o- (40) 
Thus we have proved the following important theorem. 
Theorem 4.2 Let u, u, be the solutions of (2) and (17), respectively, then we have 
|u—ualla < Chllfllo,o. (41) 


5 Enhancement of the accuracy 


We propose a simple post-processing scheme to improve the accuracy of the finite element 
solution of problem (17). Here we only state the main result concerning the mathematical 
analysis. We will present it in our another work. 

Noticing that the function va € V, on M can be written as 


9 
valu = >_ vidi, 
i=1 


where v; = v(a;) and ¢;, i = 1,2,--- ,9 are the associated basis functions, then we define an 
interpolation Hz4Va € [Q2(M)]? by 


9 
TlanValm = So viti, VME Tan; 
i=l 
where p;i, i = 1,2,--- ,9 are the associated basis functions of the biquadratic element. 
Since interpolation Il2}v} € [Q2(M)]? uses the interpolation theorem, the following theorem 
is derived directly. 
Theorem 5.1 Let u, u, be the solution of (2) and (17), respectively, then we have 


lu — zruz l|a < Ch? flao. 
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6 Numerical test 
We consider the equation (1) with w=1, f= (ft, f2) and Q = [0, 1} x [0, 1}, where 
fi = n? [4sin 2ry(—1 + 2cos 2x) — cos T(x +y) + 725 sin asin Ty], 
fo = n? [4sin 2rz(—1 + 2 cos 2ry) — cos n(z+y)+ is sin rz sin Ty]. 
Then it can be verified that the exact solution of problem (1) is u = (uy, u2) with 


u1 = sin 27y(—1 + 2cos2mx) + z4 sin mg sin ry, 


ug = sin 2rg(—1 + 2cos2my) + zh sin we sin ry. 
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The mesh on 22 is obtained by dividing it into n x n squares. As comparisons, we computed 


the energy norm of QB element and the conforming Q; element, together with |]u — Iz, ual, 


of QB element. Numerical results are shown in F igure 2 and Figure 3. It is easy to see the 


robust optimal convergence and superconvergence in the energy norm and the advantage of the 


QB element over the conforming Qı element. 
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Remark 6.1 The results in this paper can be extend to three dimensional cases with proper 


modifications. 
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